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^^ Abstract 
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^ ! Approximating complex curves with simple parametric curves is widely used in CAGD, 
CG, and CNC. This paper presents an algorithm to compute a certified approximation to 
a given parametric space curve with cubic B-spline curves. By certified, we mean that the 
approximation can approximate the given curve to any given precision and preserve the 
geometric features of the given curve such as the topology, singular points, etc. The approx- 
imated curve is divided into segments called quasi-cubic Bezier curve segments which have 
properties similar to a cubic rational Bezier curve. And the approximate curve is naturally 
constructed as the associated cubic rational Bezier curve of the control tetrahedron of a 
Q ■ quasi-cubic curve. A novel optimization method is proposed to select proper weights in the 
cubic rational Bezier curve to approximate the given curve. The error of the approximation 
is controlled by the size of its tetrahedron, which converges to zero by subdividing the curve 
segments. As an application, approximate implicit equations of the approximated curves 
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t^^ ' can be computed. Experiments show that the method can approximate space curves of high 
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degrees with high precision and very few cubic Bezier curve segments. 
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Bezier curve, cubic B-spline curve. 



^ ! 1. Introduction 



Parametric curves are widely used in different fields such as computer aided geometric 
design (CAGD), computer graphics (CG), computed numerical control (CNC) systems [1, 2]. 
One basic problem in the study of parametric curves is to approximate the curve with 
lower degree curve segments. For a given digital curve, there exist methods to find such 
approximate curves efficiently [3, 4, 5, 6]. If the curve is given by explicit expressions, 
either parametric or implicit, these methods are still usable. However, some important 
geometric features such as singular points cannot be preserved. In this paper, we will focus 
on computing approximate curves which can approximate the given curve to any precision 
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and preserve the topology and certain geometric features of the given space curve. Such 
an approximate curve is called a certified approximation. Here, the geometric features 
include cusps, self-intersected points, inflection points, torsion vanishing points, as well as 
the segmenting points and the left (right) Frenet frames of these points. 

There are lots of papers tried to approximate a smooth parametric curve segment [1, 
7, 8, 9, 10, 11, 12, 13, 14]. Among them. Geometric Hermite Interpolation (GHI) is a 
typical method for the curve approximation. Degen [8] presented an overview over the 
developments of geometric Hermite approximation theory for planar curves. Several 2D 
interpolation schemes to produce curves close to circles were proposed in [9]. The certified 
approximation were considered by some authors and they focused on the case of planar 
curves [15, 16, 17, 18]. 

For space curves, Hijllig and Koch [10] improved the standard cubic Hermite interpolation 
with approximation order five by interpolating a third point. Xu and Shi [11] considered 
the GHI for space curves by parametric quartic Bezier curve. Pelosi et al. [12] discussed the 
problem of Hermite interpolation by using PH cubic segments. Chen et al. [14] enhanced 
the GHI by adding an inner tangent point and the approximation was then more accurate. 
These methods were mainly designed for the local approximation of a parametric curve 
segment. The approximate curves obtained generally cannot preserve geometric features 
and topologies for the global approximation. The algorithms had to be improved to meet 
certain special conditions. For instance, Wu et al [19] presented an algorithm to preserve the 
topology of voxelisation and Chen et al [20] gave the formula of the intersection curve of two 
ruled surfaces by the bracket method. As a further development for certified approximation, 
more properties such as the topology and singularities of the curve need to be discussed 
in the approximation process. We would like to give the local approximation with certain 
restrictions. And the local approximation methods can then be used in the global certified 
approximation naturally. 

The certified approximation is also based on the topology determination. For implicit 
curves, the problem of topology determination was studied in some papers such as [21, 22, 
23, 24]. Efficient algorithms were proposed in [25] and [26] to compute the real singular 
points of a rational parametric space curve by the /i-basis method and the generalized D- 
resultant method respectively. An algorithm was proposed to compute the topology for a 
rational parametric space curve [27]. However, even we have the methods to determine the 
topology of space curves and the methods to approximate the space curves with free form 
curves, the combination of them is not straightforward. The topology may change while 
the line edges in topology graph are replaced by the approximate free form curve segments. 
For example, some knots may be brought in or lost such that the crossing number of the 
approximate curve is not equivalent to the approximated curve. 

In this paper, we compute a certified approximation to a given parametric space curve 
with a rational cubic B-spline curve based on the topology. The cubic rational Bezier curve 
is taken as the approximate curve segment because it is the simplest non-planar curve and 
has nice properties [28, 29]. The presented method consists of two major steps. 

In the first step, the given space curve segment is divided into sub-segments which have 
similar properties to a cubic rational Bezier curve. Such curve segments are called quasi- 
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cubic Bezier curves. The preliminary work of our division procedure is to compute the 
singular points and the topology graph of the given curve, which have already been studied 
in [30, 25, 26, 27]. Inflection points and torsion vanishing points of the curve are also added as 
character points. We further divide the curve segments to ensure that the subdivided curve 
segments have similar properties to a cubic Bezier curve. For instance, each curve segment 
has an associated control tetrahedron whose four vertices consist of the two endpoints of 
the curve segment and the two intersection points of the tangent lines and the osculating 
planes at the different endpoints respectively. And the curve segment is inside its associated 
control tetrahedron. Furthermore, we need to ensure some monotone properties about the 
associated control tetrahedron, which are necessary for the convergence of the algorithm. 
The tetrahedrons are then just the control polytope of the approximate cubic Bezier curves. 
In other words, the approximate curve is controlled by the sequence of the tetrahedrons. 
And this property ensure the topological isotopy for the approximated and approximate 
curves. Some more careful discussions are proposed for both cubic Bezier and quasi-cubic 
curve segments. 

In the second step of the algorithm, we use a cubic rational Bezier spline to approximate 
a quasi-cubic Bezier curve obtained in the first step. Some different approximation methods 
can be used here such as GHI with inner tangent points [14]. However, as we mentioned, a 
quasi-cubic Bezier curve has an associated control tetrahedron. The associated cubic rational 
Bezier curve of this tetrahedron is naturally used as the approximate curve. So, each curve 
segment and its approximated cubic curve segment share the same control tetrahedron. 
A novel method, called shoulder point approximation, is proposed to select parameters in 
the cubic Bezier curve so that it can optimally approximate the given curve segment. If 
the distance between the two curve segments is larger than the given precision, we further 
subdivide the given curve segment and approximate each sub-segment similarly. The error 
of the approximation is controlled by the size of the associated tetrahedrons, which are 
proved to converge to zero. In the subdivision process, there is one important difference 
between our algorithm with the others. We only need to check the collision of the sub- 
tetrahedrons subdivided from which are the intersected before the subdivision, since the 
sub-tetrahedrons are included in its father tetrahedrons. In general algorithms, one has to 
check the collision of all pair of the approximate curve segments or their control polytopes 
after a subdivision. Finally, the rational cubic Bezier curves are converted to a C^ rational 
B-spline with a proper knot selection and used as the final approximate curve. After a cubic 
parametric approximate segment is computed, we can compute its algebraic variety using 
the /i-basis method [31], which can be used as the approximate implicit equations for the 
given parametric curve. 

The proposed method is implemented and experimental results show that the method 
can be used to compute certified approximate curves to high degree space curves efficiently. 
The computed rational B-spline has very few pieces and can approximate the given curves 
with high precision. 

The rest of this paper is organized as follows. In Section 2, some notations and prelimi- 
nary results are given. In Section 3, we give the algorithm to compute the dividing points 
such that each divided segment is a quasi-cubic curve. In Section 4, the method of parameter 
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selection for the cubic rational Bezier segments is proposed and then an algorithm based 
on shoulder point approximation is given. We also prove that the termination of the algo- 
rithm. The final algorithm is given in Section 5, and some examples are used to illustrate 
the algorithm. In section 6, the paper is concluded. 

2. Preliminaries 

Basic notations and preliminary results about rational parametric curves and cubic Bezier 
curves are presented in this section. 

2.1. Basic notations 

A parametric space curve is defined as 

r{t) = ixit),yit),zit)), (2.1) 

where x{t),y(t), z(t) G Q(t) and Q is the field of rational numbers. In the univariate case, 
Liiroth's theorem provides a proper reparametrization algorithm and some improved algo- 
rithms which can also be found such as [30]. So we assume that (2.1) is a proper parametric 
curve in an interval [0, 1] since any interval [a, b] can be transformed to [0, 1] by a parametric 
transformation t ■^ |^. Further, the denominators of (2.1) are assumed to have no real 
roots in [0, 1]. 

The tangent vector of r(t) is r'(t) = {x' (t) , y' (t) , z' (t)) and the tangent line of r(t) at a 
point r(to) is T(to) = r(io) + Ar'(to), A G Q. A point r(to) is called a singular point if it 
corresponds to more than one parameters with multiplicities counted. A singular point is 
called a cusp if r'(to) is the vector of zeros, which means that to is a multiple parameter; 
otherwise, it is an ordinary singular point [26]. The curvature and torsion of the curve are 



|r'(t)P ' ||r' X r"|| 

A point is called an inflection if its curvature is zero and called torsion vanishing point if 
its torsion is zero. All these points are called character points of the curve, and r(t) is a 
normal curve if it has a finite number of character points. A rational space curve is always 
a normal curve. In this paper, we assume that K,{t) ^ and T(t) ^ 0, which means that the 
curve is not a planar curve. 

If r(to) is not a character point, then the Frenet frame at r(to) can be defined as J^(to) := 

{r(to); a(to),/3(to),7(to)} where a(to) = ^, /3(to) = 7(^0) x a(to), 7(^0) = \\Ct\^Sm\ 
are the unit tangent vector, unit principal normal vector, and unit bi-normal vector, respec- 
tively. And the osculating plane is 0(to) '■= {{x, y, z) — r(to)) ■ 7(^0) = 0. 

For a point with fi;(to) = 0, the bi-normal vector is not defined, neither is the osculating 
plane. Here, we define them using limit. Consider the limit limt^to 7(^) of the bi-normal 
vector at to- Since the left limit and the right limit are generally different, we define the lefl 
bi-normal vector and the right bi-normal vector as 7~(to) := limt^to_o7(t) and 7^ (to) : = 
liTat-^to+o "f (t) respectively. The limitations always exist if r(t) is a rational space curve 
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of form (2.1). As a consequence, the left and right osculating planes at to are O (to) := 
((x, y, z) — r(to)) ■ 7^ = and O'^ito) := ((x, y, z) — r(to)) ■ 7"^ = 0. If the fi:(to) 7^ 0, one can 
find that 7+(to) = 7"(to) and 0+(to) = O^(to). 

Similarly, if to is at a cusp, we define the left and right tangent vectors as a" (to) : = 
\im.t^to-o ct(t) and Q:"'"(to) := lirat^to+o ^(t) , respectively. Hence, the corresponding left and 
right principal vectors are f3~(tQ) := 7"(to) x Q;~(to) and /3^(to) := 7'''(to) x Q:+(to). We 
also denote the left and right tangent lines as T~(to) = r(to) + AQ;~(to) and T+(to) = 
r(to) + AQ:+(to) where A is the real number parameter. Then, a rational parametric curve 
r(t) always has left and right Frenet frames. 

2.2. Rational cubic Bezier curve 

A rational Bezier curve with degree n has the following form 

where Wj > are associated weights of the control points p,, G M^ and B'^if) = (") (1 — t)"'~*t*. 
When n = 3, it defines a cubic rational Bezier curve where OP0P1P2P3 is called the control 
tetrahedron of p(t). One can set the weight Wq = Wa = 1 up to a parametric transformation. 
We now consider the cubic curve and omit superscript 3 from Bf{t) 

_ poBojt) + a;ipi^i(t) + uj2P2B2{t) + psBsjt) , ^ r„ .. .. .^ 

^^' ~ Bo{t) + u,B,{t) + u^B^it) + Bs{t) ' * ^ ^'' 'J- ^'■'> 

The rational cubic Bezier curve (2.2) has the following properties. 
Lemma 2.1. Let p(t) be a non-planar cubic rational curve of the form (2.2). Then 

1) p(t) passes through the endpoints po, Ps with the corresponding tangent directions p'(0) 

and p'(l) parallel to poPi and P2P3 respectively. 

2) P0P1P2 and P1P2PZ are the osculating planes ofp(t) at the endpoints po andps, respec- 

tively. 

3) p(t) lies inside its control tetrahedron OpoPiP2P3- 

4) p(t) has no singular points and K(t) 7^ 0,r(t) ^ in [0, 1]. 

5) For any t^ < t2 G [0, 1], the control tetrahedron of p*{t) = p(t),t G [t^,t2] is inside the 

control tetrahedron of p{t) . 

6) IIpoPoiII) IIP1P12II; and IIP2P23II are strictly monotone for t* G (0,1) where Poi,Pi2; and 

P23 are the intersection points of the osculating plane 0{t*) with poPi, PiP2; and P2P3 
respectively. 

7) IIPoPosII and ||piPi2|| are strictly monotone for t* G (0, 1) where po3 = PiP2P(^*) flPoPs 

and P12 = PoP3P(^*) n PiP2- 
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Proof. Properties 1), 2) and 3) are basic properties of Bezier curves and the proof can be 
founded in [1]. They also can be checked directly. 

For 4), Li and Cripps shown that there is no cusps and inflection points for a non- 
degenerate rational cubic space curves in [32], and the torsion can be checked directly. 
Wang et al. also proved that a cubic space curve has no singular points by moving planes 
method in [25]. 

5) can be proved by a successive Decasteljau subdivision [1]. The control tetrahedron of 
Pi{t),t e [ti, 1] is inside the control tetrahedron of p(t). Successively, the control tetrahedron 
of p*(t), t G [ti, ^2] lies in the control tetrahedron of Pi(t). 

Property 6) can be derived from the above five properties. Also this property is a special 
case of the following Theorem 3.10 in this paper. 

For 7), it is sufficient to prove that the planes PiP2P(^*) and PoP3P(^*) do not touch 
p(t*) with t* G (0,1), respectively. Since PoP3P(^*) passes through pcPs and p(t) is cu- 
bic, PoP3P(^*) cannot have any tangent point different from po,P3- Supposing the plane 
PiP2p(t*) touches p(t*) at t* G (0, 1), the osculating plane 0(t*) must intersects PiP2P(^*) 
with the tangent line T(t*). By 6), T(t*) must intersect P1P2 which is the intersection line 
of 0(0), 0(1). However, according to Decasteljau subdivision, the intersection point of T(t*) 
and 0(0) is always different from that of T{t*) and 0(1). Then there is a contradiction. D 

The shoulder point of a cubic Bezier curve will play an important role [28]. The definition 
is given below. 

Definition 2.2. Let p(t) be a curve of the form (2.2). Its shoulder point s is defined as 
intersection point of p{t) and the plane P1P2PM where Pm = (Po + P3)/2 (Figure 1). 




Figure 1: Shoulder point of a Bezier cubic curve 



Proposition 2.3. Let s he the shoulder point of p{t). Then s = p(l/2) = AiPi + A2P2 
(1 - Al - A2)pm where Ai = ,,.l^'^ . As 



3u}2 



2+3wi+3w2' ^ 2+3LJi+3tJ2' 

Proof. By 7) of Lemma 2.1, there exists a unique intersection point of p(t) and the plane 
PiP2Pm- And Ai, A2 and 1 — Ai — A2 are just the area coordinates of s in the triangle PiP2Pm- 
More details can be found in [28]. D 
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It is known that the curve is closer to the control point when its associated weight 
is greater. We now consider the point which has the maximum distance to the planes 
-Pi = P0P2P3 and P2 = P0P1P3 respectively. 

Definition 2.4. Let r(t),t G [0, 1] be a curve segment on the same side of a plane Q with 
the two endpoints on Q. For another plane R parallel to Q, a tangent point ofr(t) with the 
plane R is called a parallel point ofr(t) associated to the plane Q. 

According to the definition, a parallel point should satisfy 

|r'(t),qi -qo,q2-qo| =0, (2.3) 

where qo, qi and q2 are three non co-linear points on Q. In general, there may be several 
parallel points for a curve segment and a fixed plane. However, for the rational cubic curve 
segment (2.2), there is a unique parallel point associated to Pi = P0P2P3, and similarly, 
there is a unique parallel point associated to P2 = PoPiPs- 

Proposition 2.5. Let p(t) be a non-planar cubic rational curve of the form (2.2). Then 
there are unique parallel points associated to the planes Pi and P2 respectively, and they are 
points of p{t) having the maximal distance to Pi and P2 respectively. 

Proof. By equation (2.3), we can find that ^^^~^^*^+ff~^ = Ui and ^^^ ^tfutl)"^ = ^2 are the 
constraint equations for the parallel points associated to Pi and P2 respectively. They are 
two monotone functions for t E (0, 1) with two asymptotes t = 0, 1. It means that for any 
weights there is only one parallel point associated to Pj. Furthermore, the parallel point has 
the maximal distance since the endpoints of the curve are on Pj. D 

3. Quasi-cubic segments on space parametric curves 

In this section, we propose a method to divide a given curve r(t) into segments which 
have similar properties to cubic Bezier curves, which are called quasi-cubic Bezier segments 
and can be approximated by cubic rational Bezier curves nicely. 

3.1. Conditions for subdivision 

Let to, ti be the endpoints of a curve segment r(t). We will define an associated tetrahe- 
dron for it. Let 0~^(to) and 0~(ti) be the right and left osculating planes at the endpoints 
respectively. We denote their intersection line as L, if they are not parallel. Since L and the 
right tangent line T''~(to) are coplanar, they intersect at a point ri if they are not parallel. 
Similarly, L and the right tangent line T~(ti) intersect at a point r2 if they are not parallel. 
So we obtain an associated tetrahedron (}{to,ti) = Oi"orir2i"3 where ro = r(to) and r^ = r(ti) 
if ri ^ V2. 

We have shown that a cubic Bezier curve segment has eight properties in Lemma 2.1 
and Proposition 2.5. In the following, we will show how to divide any given rational curve 
segment into sub-segments having similar properties. 
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Definition 3.1. A curve segment is called a quasi-cubic Bezier curve segment, or simply a 
quasi-cubic segment, if it has the eight properties in Lemma 2.1 and Proposition 2.5. 

Theorem 3.2. Given r(t) and to, there always exists ti > to such that r(t),t G [to,ti] is a 
quasi-cubic Bezier curve segment. 

We leave tlie proof of tliis tlieorem at tlie end of tlie subsection 3.2. 

Definition 3.3. Let r{t),t G [toi^i] be a quasi-cubic segment. Then its associated cubic 
Bezier curve segment is defined by the associated tetrahedron ofr{t), i.e., the control points 
are Fq, ri, r2 and r-^. 

In order to divide the curve segment into quasi-cubic segments, we first add the inflection 
points and torsion vanishing points as the dividing points, denoted by P. The parameters 
of these points can be computed by solving the real roots of K(t)r(t) = 0. The left and right 
Frenet frames are also needed. There are several efficient methods to find the real roots of 
a univariate polynomial [33, 34] and one can use the procedures realroot and isolate in 
Maple. 

We need to find more dividing points. Fix a start point t = to, we now try to determine ti 
such that ti—to is as big as possible and the segment is included in its associated tetrahedron 
designed above. Several boundary parametric values to exclude some special points with 
respect to to are computed in the following cases: 

Condition I). Let t\ > to be its nearest parametric value from P. Find ti G (to,t*) 
such that Fi{si,S2) := Q:"''(si) ■ 7~(s2) 7^ and F2{si,S2) := Q:~(s2) ■ 7'^(si) 7^ for any 
to < si < S2 < ti, meaning that the right tangent vector a"'"(si) is not parallel to the 
left osculating plane 0~(s2) and the left tangent vector a~(s2) is not parallel to the left 
osculating plane 0"'~(si). 

Since the curve is non-planar, Fi{si,S2),i = 1,2 cannot be identically zero. We take 
a further look at the inequalities Fi 7^ 0, F2 7^ 0. Since the derivative can be computed 
using limits, r(t) is differentiable to any order although the left and right derivative may be 
different. For conveniences, we omit the -|-, — marks to distinguish between left and right 
derivatives. In what below, we give detailed analysis for Fi and the analysis of F2 is similar. 

Fi{si, S2) = cx{si) ■ 7(52) = II ,, .|||| ,, . ,,, .|| . 

||r'(si)||||r'(s2) X r"(s2)|| 

Assuming Si = to + Si, S2 = Si -\- 62, Si > 0,62 > 0, Fi{si, 82) is re-parameterized as 

|r'(to + 61), r'(to + Si + 62), r"(to + Si + 62)1 



Fi{6i,62] 



|r'(to + 5i)||||r'(to + 5i + 82) x r"(to + 5i + 82) 



Expanding the vectors of the numerator at t = to+Si as Taylor series r'g(to+Si), r'g(to+Si-\-62) 
and r'^(to -{-61 + 62) respectively, and combining them, we have 

p.. _ 6l\r'sito + ^1), F^(to + ^1 + ^2), r'lJito + 61 + ^2)! , . 

'^ '' '' \\r'{to + 6i)\\\\r'{to + 6i + 62)xr"{to + 6i + 62)r ^ ' ^ 
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where f'^(to + ^i + ^2) = (r^(to + ^1 + ^2) - r'gito + 6i))/S2 and f'^'(to + ^i + ^2) = (r"(to + 
^1 + ^2) - rg{to + Si + S2))/S2. Furthermore, when 62 = 0, r'gito + 61) = r'sito + 61) and 
^'"(t + S ) = r"'(t + S ) 
' LetM6uS2) = F,{S,,62)/6l Then f,{6u0) = r{to + 6,)/\\r'{to + 6^)\\. Fi(5i, ^2) = is 
a planar curve in the plane of (5i, ^2) which has two components: a double line 5| = and 
another planar curve f 1(61,62) = 0. That means /i = intersects ^2 = with the points 
which are exactly the torsion vanishing points r(to + 5i) = of r(t). And we need not 
compute these points since they are already included in the separating points needed in the 
topology computation which is discussed in Section 3.3. Consider the intersection points of 
fi{6i,62) and Si = 0. We can find that the real roots of fi{0,62) = are associated to the 
vector a(si) = r'(to) just parallelling to the osculating plane 0(^2) = 0(to + 62)- 
Thus, condition I) can be reduced to solve the following optimization problem 



min 61 + 62 

s.t. Fi{5i,S2) = 0,5i>0,S2>0 



(3.2) 



and then ti can be selected from (to, to + 61 + 62). There are numerical methods to solve the 
optimization problem. However, we prefer to solve it based on the above discussion since 
it is enough to get a boundary parametric value less than the exact solution of (3.2). We 
can find the positive real roots of /i(5i, 0) and /i(0, ^2) for 61 and 62 respectively. Let 6^ be 
the minimal one among all the real roots. Then 61 + 62 = S* defines a line. If the line does 
not intersect /i in the first quadrant, then ti can be in (to, to + S*). This can be checked by 
finding the real roots of fi{Sl — 62, 62) = 0. Otherwise, set S* -^ S*/2 and check the process 
repeatedly until the proper 6* is found. If /i(5i, 0) and /i(0, ^2) have no positive real roots, 
6^ can be initialed as 6* = t* — tQ. 

Similarly, we can find such a ^2 for F2. Finally, let t2 = min(to + 5i,to + ^2) t>e the 
boundary parametric value of ti. 

Remark 3.4. The function Fi{6i,62) in (3.1) actually has a finite number of terms if the 
approximated curve r is a rational curve. If r is a parametric curve in elementary functions, 
-^1(^17^2) will be in the series form. However, the problem (3.2) can still be solved using 
a numerical method. Starting with an initial value 6^, we can find a boundary number by 
checking whether 61 + 62 = SI and Fi{6i, 62) have common points in the first quadrant with 
one of the directions {5^ -^ '^i/2, 6^ •^ 2(5^}. 

Further restrictions will be proposed afterward. We will omit the similar discussions and 
solving processes and give the conditions directly. 

Condition II). Let tg be the parametric value ti computed in the above procedure. 
Find ti G (to, t2) such that 

F{si, S2) := a+(si) x (r(s2) - r(si)) ■ a-(s2) ^ 

for any to < Si < S2 < ti, which means that the right tangent line T"''(si) and the left 
tangent line T~(s2) are not coplanar. 
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Condition III). Let ^3 be the parametric value ti computed in the above procedure. 
We should find ti G (^0,^3) such that -Fi(si,S2) := 0^{s2){r{si)) 7^ and ^2(51,52) := 
0"*"(si)(r(s2)) 7^ 0, which imply that r(si) is not on the left osculating plane 0~{s2) and 
r(s2) is not on the right osculating plane 0~^{si). 

Conditions I), II), and III) are used to guarantee that the tetrahedron Oi"oi"ii"2i"3 is not 
degenerated to a plane polygon. However, these conditions are still not sufficient for the 
curve segment lying inside Ororir2r3. We will give one more condition such that the curve 
segment lies inside the tetrahedron and has only one parallel points associated to planes Pi 
and P2 respectively. 

Let ti < ^4 where ^4 is the parameter value obtained from III). Then the curve segment 
r{t),t G [to;^i] satisfies the conditions of I) to III) and r(t) has no character points. We will 
try to find t* G (to,ti] such that for any Si < S2 < Sa G [to,t*], the tangent vectors a(si), 
Q:(s2), and 01(53) are not coplanar, i.e., 

\ol{si),oc{s2),oc{s3)\j^0. (3.3) 

The following lemma is needed for further discussion. 

Lemma 3.5. For a fixed to and Ve > 0, F{si, S2) '■= \(y.{to), Q:(si), 01(52) | = has solutions 
(si, S2) in (0, e)^ if and only if r(t) is a planar curve. 

Proof. It can be checked by expanding vectors to Taylor series which are partly illustrated 
above. D 

And the lemma also holds for F mentioned in I) to III). It means that F{si,S2) has no 
branch segment on the first quadrant of the (si, S2) plane connecting the origin point. 

Condition IV). Find t* G (to,ti) such that F := |q;(si), 0(32), 0(53)! ¥" fo^ a,ny 
Si < S2 < S3 G [to,t*] C [to,ti]. That means r(t) does not have a triple of linear dependent 
tangents in [to, t*]. Suppose si = to + 61, S2 = si + S2 and S3 = S2 + S3 where 61 > 0,62 > 
and ^3 > 0. 

If ^1 > 0, then we need to find the least to + Si + 62 + S3 with F{Si, S2, S3) = 0, that is, 

min Si + S2 + S3 
s.t. F{Si,S2,S3) = 0,Si,S2.S3>0. 

By Taylor expansion, we find that F((5i, ^2, ^3) has no branch passing through the {Si,Sj) 
plane from the first octant in the space of {Si, S2, S3). Then we initialize Si, i = 1, 2, 3 in the 
plane S1 + S2 + S3 = Si = ti and check the intersection of the plane with F. Set the boundary 
parametric value tg^ = Si if there is no intersection; otherwise set Si ■^ S\/2 and repeat the 
checking process. 

If Si = 0, then F{S2, S3) degenerates to the special case mentioned in Lemma 3.5 and we 
can find a boundary parametric value as ^52- Finally, let t* = iniii{tt,i,t'^2)- 

We have the following key theorem. 

Theorem 3.6. Let t* he found by the above process. For any e > 0, ti = t* — e > to, the 
associated tetrahedron (>i"orir2r3 ofv{t),t G [to7^i] ^s not degenerated. Furthermore, 
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1) r(t) passes through the endpoints ro,r3 with the corresponding tangent directions r'(to) 

and r'(ti) parallel to rori and r2r3 respectively. 

2) rorir2 and rir2r3 are the osculating planes ofr{t) at the endpoints tq andr^, respectively. 

3) r(t) lies inside its control tetrahedron Oi"orir2r3. 

4) r(t) has no singular points and nit) ^ 0,r(t) ^ in [to,^i]- 

5) There exists only one parallel point between ri and ror2r3, same to r2 and rofirs. 

Proof. According to conditions I) to III), the tetrahedron <0i"orir2r3 does not degenerate. 
1), 2), and 4) are also followed by the discussions. 

The curve segment is inside the tetrahedron. We claim that the curve segment and r^ are 
on the same side of plane P3 = rorir2. Otherwise, there exists a parallel point p associated 
to P3 but on the different side with r^, since r(t) is a smooth segment. Then cx.{p) is parallel 
to P3 which contradicts to I). Similarly, the curve and fq are on the same side of Pq = rir2r3. 
Furthermore, the curve and ri are on the same side of Pi = ror2r3. Otherwise, there exist at 
least two parallel points pi, p2 on different sides of Pi. Then |q;(pi), Q:(p2), Q:(r3)| = which 
contradicts to condition IV). Similarly, the curve and r2 are on the same side of P2 = rorir3. 
Therefore, 3) is followed. 

Finally, 5) is correct. Otherwise, there exist at least two parallel points associated to Pi 
or P2 which will lead a contradiction to condition IV). D 

Proposition 3.7. For any t^ < tg ^ [^Oj^i]; the sub-tetrahedron Oro''i''2'^3 ^/ the sub- 
segment r*(t),t G [^1,^2] also has the properties listed in Theorem 3.6. 

Proof. In the dividing process, the conditions in I) to IV) are satisfied for the parameters 
through the interval not just only for the endpoints. Then the properties are all satisfied 
within [t\,t*^] C [to,ti]- □ 

3.2. Further properties of the divided segment 

In this subsection, we prove that the curve segment obtained in the preceding section 
also has properties 6) and 7) in Lemma 2.1. Before that, we need some preparations. 

Suppose that the curve segment r(t),t e [^07^1] satisfies conditions I) - IV) in the pre- 
ceding section. 

Lemma 3.8. Let Oi"oi"ii"2r3 be the control tetrahedron of a given curve segment r(t),t G 
[to,^i]- Then for any t* G (to,ti), the control tetrahedron (>i'orii"2''3 ^f the curve segment 
r{t),t G [to,t*] has the following properties: 

1. r\ and ri are on the same side of tq in the tangent line T(to); 

2. r2 and r2 are on the same side o/T(to) in the osculating plane 0(to)- 

Proof. Using the first and second order Taylor expansion of r(t), one can prove the lemma. 

D 
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Lemma 3.9. Let 0(t*) be the osculating plane of curve rit) at t* G [to,ti]. Ifr(t) does not 
pass through 0{t*), then T{t*) = 0. 

Proof. Similar to the discussions of condition I), using the third order Taylor expansion, 
one can see that |r'(t*), r"(t*), r"'(t*)| = 0, that is r(t*) = 0. D 

We now prove another key property for the curve segments. 

Theorem 3.10. Let Oi"orir2i"3 be the associated tetrahedron of a curve segment r(t),t G 
[to,ti]. Then ||roroi||, ||riri2||, and ||r2r23|| are strictly monotone in {tQ,ti) where roi,ri2, 
and r23 are the intersection points of the osculating plane 0{t*) and rori,rir2, and r2r3 
respectively. 

Proof. Firstly, the intersection point roi of rori and the osculating plane 0{t*) must be 
on the same side with ri with respect to ro on the curve segment. Otherwise, subdividing 
r(t) at t*, the sub-segment r^(t),t G [to,t*] will not be inside its tetrahedron for fqi 7^ fq by 
Lemma 3.8. We denote by ro2 the intersection point of line ror2 and 0{t*). Similarly, r23 
is on the same side with r2 with respect to rs and ro2 is on the same side with r2 w.r.t. tq 
(See Figure 2). 



Figure 2: The osculating plane 

Secondly, we claim that there exist no t\ < tg in [^o^^i] such that the osculating planes 
0(t\) and 0(^2) have the same intersection point roi with rori. It is sufficient to prove that 
there has no t* G (to, ti) such that the osculating plane 0(t*) passes through ri by assuming 
^2 = ti and denote t* by t*. Otherwise, if the osculating plane 0{t*) passes through ri, then 
0(t*) passes through the line rir(t*) but cannot pass through fq and r3 by the restrictions 
in condition I). Hence 0{t*) has only two possible cases: it either intersects rora and the 
polygonal line ror2r3, or intersects ror2 and r2r3. In the first case, let the intersection points 
of T(t*) and 0(to), 0(ti) be foq, roi respectively. Then roo and tq^ are on the same side with 
respect to r(t*) in line T(t*). Which means that one of the sub-segments r^(t),t G [to?^*] 
and r2(t),t G [t*,ti] cannot be inside its tetrahedron by the first paragraph of the proof, a 
contradiction to Proposition 3.7. In the second case, the points ro and r3 are on the same 
side of 0(t*). By Proposition 3.7, the sub-segment curves at t = t* are also on the same side 
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of 0{t*). Then the curve r(t) does not pass through 0(t*) at t*, which means that r(t*) = 
by Lemma 3.9. Hence, ||roroi|| and ||r2r23|| are monotone. 

It is known that Fqi hes on rori and r23 hes on r2r3. We claim that ri2 must be on rir2. 
Otherwise, assuming 0{t*) has no common points with rir2, then 0(t*) must intersect with 
rori,ror2, FiFa, and r2r3. That means fq and r3 are on the same side of 0(t*), and then 
r(t*) = 0, a contradiction. 

Since the curve is inside its tetrahedron, r(t*) is inside the quadrangle roiri2r23r3o. Ac- 
tually, r(t*) is inside the triangle roiri2r23. r(t*) cannot be on roi and r23 according to 
condition III). So, if r(t*) is not inside the triangle roiri2r23, then r{t*) is on the opposite 
side with ri2 with respect to roir23 or on roir23. Then T(t*) f]0(to) is not inside roiri2, or, 
T{t*) f]0{ti) is not inside ri2r23, since roiri2r23r3o is convex. Without loss of generality, 
we suppose T{t*)f]0{to) is not in roiri2. Then T(t*)f]0{to) are on the same side with 
r2 w.r.t. rori in 0(to) by Lemma 3.8. Hence, T{t*)f]0{to) and T{t*)f]0{ti) is on the 
same side of r{t*) in T(t*), which means that one of the sub-segments r^(t),t G [to;^*] and 
r2(t),t G [t*,ti] cannot be inside its tetrahedron, a contradiction to Proposition 3.7. 

Therefore, r{t*) can only be inside the triangle roiri2r23, and T(t*) can only intersect 
i"oii"i2 with roi2 and intersect ri2r23 with ri23. Subdivide r(t) at t = t* to get curve seg- 
ments r^(t),t G [to,t*], and r2(t),t G [t*,ti], and their tetrahedrons as C>roi'oiroi2i"(^*) and 
(>r(t*)ri23r23r3. It has been shown that these two sub-tetrahedrons are inside the tetrahe- 
dron Oroi"ir2i'3- As a consequence, for any t^ < ^2 in [to^^i], the sub-tetrahedron of the 
sub-segment r*{t),t G [^^,^2] is inside the tetrahedron <>i"orir2r3. 

Finally, we prove that ||riri2|| is monotone. It suffices to show that there exist no 
t* < tg G [to,ti] such that O(t^) and 0(^2) have a common point in rir2. Otherwise, 
we assume O(t^) and 0(t2) have a common point r^2 ^^ rir2. Since roFoi and r2r23 are 
monotonously increasing, roi(ti) and r23(t^) are on the same side of 0(^2)- Hence the 
intersection line of O(t^) and 0(^2) can only be outside of the tetrahedron <)roi"ii"2i"3 passing 
through r^2- Then the sub-tetrahedron of the sub-segment r\2it),t G [t\, tg], cannot be inside 
the tetrahedron Ororir2r3, which contradicts to the consequence in the preceding paragraph. 

D 

For clarity, we summarize the properties mentioned in the proof of the above theorem as 
follows. 

Proposition 3.11. For any t\ < t2 ^ [to,ti\, the sub-tetrahedron Orgr^rgrg of the sub- 
segment r*{t),t G [^1,^2] is inside the tetrahedron <)i"oi"ir2r3. 

Similar to 7) of Lemma 2.1, we have the following proposition. The proof is also similar 
to that of 7) of Lemma 2.1. 

Proposition 3.12. ||roro3|| and ||riri2|| are strictly monotone with t* G {to,ti) where ro3 
and ri2 are the intersection points rir2r(t*) f]ror3 and ror3r(t*) f]rir2 respectively. 

Proof. It is sufficient to prove that the planes rir2r(t*) and ror3r(t*) are not tangent to 
r(t) at t* G (to, ti). If the plane rir2r(t*) is tangent to r(t) at t* G {to, ti), then the osculating 
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plane 0(t*) must intersect rir2r(t*) with the tangent hne T(t*). By Theorem 3.10, T(t*) 
must intersect rir2 which is the common hne of 0(to) and 0(ti). Dividing the curve segment 
into two sub-segments r^(t) and r2(t), then one of them cannot be inside its sub-tetrahedron 
according to Lemma 3.8 which contradicts to Proposition 3.11. And one can similarly discuss 
the case for the plane ror3r(t*). D 

According to Proposition 3.12, r(t) and the plane rir2rA,/ have a unique intersection point 
Sr where Tm = (i"o + i'3)/2- We call Sr the shoulder point of the segment r(t),t G [to,ti]. 
Similar to Proposition 3.7, we can see that Theorem 3.10 and Proposition 3.12 also hold for 
any subsegment r*{t),t G [ti,t2]- 

When we subdivide the approximated curve segment at a point t = t*, by Theorem 3.10, 
we assume that the osculating plane 0{t*) intersects rori,rir2 and r2r3 at roi,ri2 and r23 
respectively. Then, one can have the following corollary. 

Corollary 3.13. Let ki{t*) = ig^, k2{t'') = ig^ and ks{t*) = i^, then ki{t*) is mono- 
tone and ki{t*) G (0, 1) with t* G (to,ti), i = 1,2,3. 

We finally give the Proof of Theorem 3.2 by summarizing the above discussions. 

Proof. Set ti as Theorem 3.6, then r(t), t G [to, ti] has the eight properties in Theorem 3.6, 
3.10 and Propositions 3.11, 3.12. It means that the segment r(t),t G [^0,^1] is a quasi-cubic 
segment. D 

3.3. Subdivision algorithm 

As we mentioned in the introduction, the topology graph ^ of a parametric space curve 
can be computed by the method in [27]. 

A topology graph is a graph Q = {V, S} where V is a set of points in the Euclidean space 
V = {vj = (aj,/3i,7j)} and £^ is a set of edges S = {(vj, Vj)|vj, Vj G V}, any two edges do 
not intersect except in the endpoints. A graph ^ is a topology graph of a parametric space 
curve r(t) if Q and r(t) have the same topology. 

The singular points of the space curve are included as vertices in Q. In this paper, we 
need to add more information to the vertices in our algorithm. For each vertex Vj in the 
topology graph, we now update it to 

Vi = {vi = r{tiQ),{tio,tii, . . . ,tik}, 

{^^v.■■■.^^k}AJ'to.■■■.J'tk}}^ (3-4) 

where each tjj is a real parameter such that r(tjj) = Vj, J-^~ and T^ are the left and right 
Frenet frames of Vj with respect to the parameters tij,j = 0, . . . ,k. The point set V thus 
updated is called the extended vertex list. Methods to compute the limitation of the tangent 
are also introduced in [23]. 

The edges in Q are not used directly in our approximation algorithm, but they give the 
connection relationship of two updated vertices. Since the space curve is parametric, the 
connection relationship is given by the parameters corresponding to the points in V in the 
increasing order. So in our paper, we use the extended vertex list V instead of topology 
graph. 
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Example 3.14. Figure 3 (a) shows a space curve with a cusp, whose topology graph is 
given in Figure 3 (h). Figure 4((i) shows a numerical approximate curve which does not pass 
through the cusp. We may use the topology graph or a refined topology graph to approximate 
the curve segment as shown in Figure 4(b)- This method has two drawbacks. First, we 
generally needs hundreds even thousands line segments to approximate the curve segment for 
a small precision [24]. Second, the approximate curve cannot keep the tangent directions of 
left and right sides of the cusp point. In this paper, we use a cubic Bezier curve instead of a 
line segment as shown in Figure 4(c); which is not only more precise but keeps the geometric 
properties of the original curve. 





(a) Origin curve (b) Topology graph 

Figure 3: Topology graph of the curve 






(a) General numerical method (b) Based on topology (e) Proposed method 

Figure 4: Numerical approxiuiate curve 

Based on the above analysis, we now give the segment dividing algorithm. 

Algorithm 3.15. Curve Subdivision. 

Input: A normal curve segment r{t),t G [0, 1]. 

Output: An extended vertex list with elements as (3.4). 
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1. Compute the certified vertex list V with all character points as vertices with the 
method in [27]. The parameters and the left and right Frenet frames are recorded. 
Suppose the real roots associated to the character points are Si,i = 1, . . . ,1 — 1 and 
= So < Si < ■ ■ ■ < s/ = 1. 

2. Divide each interval [sj, Sj+i] as Sj = Sjo < Sj^i < ■ ■ • < Sj^fc. = Sj+i such that each 
segment satisfies the conditions given in I) to IV). 

3. Rearrange the Sij in an ascending order and rename them as ti,i = 0, . . . ,n. Find the 
left and right Frenet frames of each segment r(t),t G [ti,tj+i]. 

4. Add all these new points to the extended vertex list V which is now ready for approx- 
imation. 

Each curve segment is defined by two adjoint vertices of V. By Proposition 3.7, the curve 
segment from the algorithm is in the tetrahedron and has the properties in Theorems 3.6, 
3.10 and Propositions 3.11, 3.12. Hence each curve segment obtained from Algorithm 3.15 
is a quasi-cubic segment and so are its sub-segments. 

4. Shoulder point approximation 

In this section, we propose an efficient algorithm to construct a set of cubic Bezier 
curve segments which approximate a quasi-cubic segment obtained in Algorithm 3.15 to any 
approximate bound. 

Firstly, we focus on one quasi-cubic segment r(t),t G [to,ti]. Let ro,r3 be the endpoints 
of the segment, ri the intersection point of the tangent line at ro and the osculating plane 
of r3, and r2 the intersection point of the tangent line at r3 and the osculating plane of tq. 
Then {tq, ri, r2, r3} defines a family of rational cubic curves 

/ N roEo(s) + uiriBijs) + a;2r2^2(g) + rsBsjs) 

Bo{s) + ujiBi{s} + U2B2{s) + Bs{s) 

Then p(ci;i,a;2, s) is called the associated cubic Bezier curve segment of r(t). It has been 
shown that p(ci;i, Co'2, s) meets r(t) at its endpoints r(to) and r(ti). Furthermore, p(co'i, c<;2, s) 
and r(t) have the same left and right tangent directions and osculating planes at the end- 
points, and the same control tetrahedron Oi"oi"ii"2r3. 

Proposition 4.1. Let p{uJi,uj2,s),s G [0,1] be the associated cubic Bezier curve segment 
ofr{t),t G [to5^i]- Then p{uJi,ijj2,s) can approximate r(t) at their endpoints with order two 
by setting proper ui and U2, i.e., {p(0) = r(to),p(l) = ^(ti)} and {p'(0) = r'(to),p'(l) = 
r'(ti)}. 

Proof. Following the construction of p(s) for r(t), they are G^ interpolated at their end- 
points with arbitrary weights Ui and U2- According to the properties of the cubic Bezier 
curve, one can set the proper ui and U2 such that p(s) and r(t) are C^ interpolated at their 
endpoints. D 
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In Proposition 4.1, the weights are selected to enhance the approximation order from 
G^ to C^ at the endpoints. Actually, on can get {p((X'i,a;2, 0) = r(to),p{ui,U2,l) = r(ii)} 
and {p'{ui,U2,0) = kiUir'(to),p'{ui,U2,l) = k2UJ2r' (ti)} , where ki and /c2 are positive 
constants. Hence we can set Ui and U2 such that kiUi = 1 and k20J2 = 1- However, in the 
following paragraphs, we would like to use the freedom of weights to minimize the position 
approximation error. Hence, we will show how to compute the proper weights u!i,U2 such 
that p(s) is an optimal approximation to r(t). 

The selection of the weights often leads to some optimization problems such as Tairii^-^j^j^ 
(jnaxs^t d{uji, U2, s, t)^) where d{ui, U2, s, t) is the distance function between p(a;i, a;2, s) and 
r(t) in certain forms [3]. The computation is usually not efficient and some global error 
analysis is introduced to simplify the optimization problem [35]. Another possible method 
is to approximate the target curve segment by checking the parallel points. We can push 
the parallel points of the approximated curve and the approximate curve (4.1) as near as 
possible. It also leads to an optimal problem for a function with degree three. In the 
following, we introduce a novel method which avoids any optimizations. 

The shoulder point Sp of p(s) is given in Proposition 2.3. The shoulder point Sr of r(t) 
can be computed as the unique intersection point of r(t) and the triangle rir2rj\/. Supposing 
the plane P{x, y, z) is defined by ri, r2, and tm, then the shoulder point corresponds to a real 
root t* e {to,ti) of Por(t) with r(t*) lying in the triangle rir2rA'/- So _D(a;i,a;2) = ||sp — Srp 
is a rational function in ui, U2 with total degree two. Finding the positive solution from the 
equations 

dD 



duji 
dD 



0, 

(4.2) 

0, 



we obtain the weights for the approximate cubic curve (4.1). 

Before the approximation, we will estimate the error between the two curves. Since there 
does not have any simple method to compute the distance of two parametric curves with 
different parameters, we use the distance between r and the implicit variety of a rational 
cubic curve p. It has been proved that the associated implicit ideal Jp of p can be computed 
using the /i-basis method [31] efficiently: 

Lemma 4.2. The associated ideal of p has the form Ip = {f{x,y,z),g{x,y,z), h{x,y,z)), 
where f,g and h are quadratic polynomials, i.e., the resultants ofp's fi-basis in pairs. 

The algorithm of /i-basis is given in [36]. Generalizing the approximation error function 
in [37], we have 

/ fM2 X 1/2 



,/,.(r)2 + /,(r)2 + /,(r)^ 

Let e(p, r) := e(/, r) + e{g, r) + e{h, r) = e{t) be the univariate error function in t. Then the 
approximation error can be set as the following optimization problem: 

e = max (e(t)). 

to<t<ti 
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There are many methods to solve this problem. However, for the efficiency in practice, we 
often sample t as tj = ^~ ° ,i = 0, . . . ,m, for a proper m, say m = 300, and set the 
approximate error as max(e(tj)). 

The following algorithm is proposed to approximate a quasi-cubic curve segment via 
shoulder point approximation. 

Algorithm 4.3. Shoulder point approximation 

Input: A quasi-cubic curve segment r(t),t G [to)^i] (^iT'd a positive error bound 6. 

Output: A set of cubic Bezier curves which is a 6 -approximation for r(t). 

1. Construct the associated tetrahedron of r(t) and the rational Bezier cubic curve 
p{uJi,U!2, s), s G [0, 1] as shown in (4.1). 

2. Compute the weights (cji,C(;2) such that ||Sp — Sr|| is as small as possible. 

(a) Compute shoulder points Sr and Sp(ci;i,u;2) of r(t) and p(s) respectively. 

(b) Find a pair of real roots {001,002) by solving the equation system (4.2). 

3. Compute the approximate error 6 = e(t). Ii6<6 then output p(s). Otherwise, divide 
r(t) to two parts on its middle point of arc length and repeat the approximation process 
for each subsegment. 

Example 4.4. A curve segment r(t), t G [0, 21/32] represented by the black curve with degree 
six is given by Algorithm 3.15 and the approximate cubic Bezier curve is the red dash curve 
in Figure 5. The weights are ooi = 002 = 1 in the left figure. After executing step 2 of 
Algorithm 4-3, we have ooi = 5/11, aJ2 = 16/31 in the right figure. The numerical errors are 
0.29 and 0.04 respectively computed from error function e{t) by setting m = 300. 





Figure 5: Selecting the weights for Bezier cubic curve 



To show the termination of the above algorithm, we need the following lemma. 

Lemma 4.5. The edge of the sub-tetrahedron in Algorithm 4-3 converges to zero when the 
arc length of its subdivided curve segment converges to zero. 
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Proof. There exists a t = t* G (to,ti) such that ki = 1/2 since ki{t) is monotone with t in 
(to,^i) by Corollary 3.13. Consider the subsegment r{t),t G [to,t^] and subdivide it at t = ^2 
such that k2 = 1/2 for the sub-tetrahedron 0(^0, ^i)- Then, subdivide r(t),t G [to, ^2] at t = tg 
such that ks = 1/2 for 0(to, ^3). Let t^-^^ = t^. We obtain a subsegment r(t), t G [to, t^^-*] whose 
sub-tetrahedron <)(to, f-"^^) has vertices Fq = ro,r^ ,r2 ,r3 . Similarly, we can construct 
i"j , J = 0, 1,2,3 and t^*^ According to the subdividing process, let r^- = rj,j = 0, 1,2,3. 
Then, we have H^'^rPH < \\rt'^rt'^\\/2, \\r?4^\\ < \\rt'^4^'^\\/2+\\4^'\t'^\\/2 and 
l|r?r?ll < l|rr'^rf"'^||/2 + \\r'^^~^\t'^ \\ + ||4'~'^4'"'^||/2 for i > 0. Hence, the lengthes 
of the three edges ||rQ r^ ||, ||r^ r2 || and ||r2 ''s II of ^ sub-tetrahedron 0(to,t'^*^) converge 
to zero when i -> 00. Since r(t),t G [to,ti] is a rational curve and has no singular point, 
t^*-* — to converges to zero when i —^ 00. 

Lett G [to,ti] andOror'i(t)r'2(t)r^(t) its tetrahedron. Thens(t) = ||ror;(t)|| + ||r;(t)r^(t)|| + 
||r2(t)r3(t)|| converges to zero when t — ;■ to, since r(t) has no singularities in [to,ti]. Hence 
when the arc length of its subdivided curve segment converges to zero, which means t — )■ to, 
the edge of sub-tetrahedron converges to zero. D 

The termination of Algorithm 4.3 can be guaranteed by the following theorem. 

Theorem 4.6. In Algorithm 4-3, the approximation error converges to zero for the subdi- 
vision procedure. 

Proof. By Lemma 4.5, when the arc length of its subdivided curve segment converges to 
zero, the edge of the sub-tetrahedron converges to zero. Since the approximation error is 
controlled by the edges, it converges to zero for the subdivision procedure. D 

Remark 4.7. In Algorithm 4-3, the Step 3 is given to simplify the proof of the convergence. 
In fact, for less computation, we always implement the algorithm with the following step 
instead of 3. 

3' . Compute the approximate error 6 = e(t). If6<6 then output p(s). Otherwise, divide 
r(t) to two parts on its shoulder point Sr repeat the approximation process for each 
subsegment. 

According to the proof of Lemma 4-5, the algorithm fails if a subsequence of Si does not 
converge to zero under shoulder point subdivision process, and it never happened in our 
experiments. It is an interesting problem to prove the termination of this version of the 
algorithm. 

5. Algorithms and experimental results 

After dividing the curve to segments by Algorithm 3.15, we can approximate each curve 
segment by the shoulder approximation method in Algorithm 4.3. In this section, we give 
the main approximation algorithm and the experimental results. 
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The global approximation is based on the local approximation and topology determi- 
nation in the above sections. Some relationships of the approximate curve segments are 
considerable in the global view. In our approximation, the line edges in the topology graph 
are replaced by the associated cubic Bezier curve segments. To ensure the topological iso- 
topy before and after the replacement, we restrict the cubic curve segments to have the 
appropriate topology based on the topology graph. 

It is shown that an associated cubic Bezier curve segment are decided by its tetrahedron. 
Let OP0P1P2P3 ^^^ OP0P1P2P3 be two control tetrahedrons of two cubic Bezier curve seg- 
ments p^(s) and p^(s). Then p^(s) and p^(s) can have no common points except for their 
endpoints. In the further consideration, we give two cases for the problem. The first case is 
that p^(s) and p^(s) have only one common point being the endpoint and the same Frenet 
frames at this endpoint. And the other positional situations of p^(s) and p^(s) are included 
in the second case. 

If all the pairs of cubic Bezier curves satisfy the second case, then to ensure that cubic 
curve segment does not bring in the unexpected knots while it replaces the line edge, one 
can give a sufficient condition that each cubic curve segment has no common points with 
the control tetrahedron of another curve segment except for the endpoint. This condition 
can be strengthened if we do not want to check the collision between a cubic curve segment 
and a tetrahedron. The condition can be that the two tetrahedrons have no inner points. 
By Lemma 4.5, the condition can be satisfied by subdividing the curve segments. Then the 
approximate curve have same topology with the given curve, since the approximate curve 
is controlled by the sequence of the tetrahedrons. Each tetrahedron has no common inner 
points with other tetrahedrons. 

We then only need to discuss the pairs of cubic Bezier curves belong to the first case. 
Assuming pj = Pq, then p^ is on the radial (1 — A)pj!, + Ap}, A > 0, and pg is on the same 
side with p^ on the plane pjpjpg. According to the monotonicity of the Bezier curve in 
Lemma 2.1, p^(s) and p^(s) can replace the their associated line edges without topology 
modification. 

Algorithm 5.1. Certified B-spline approximation with error bound. 

Input: A normal curve segment r{t),t G [to;^i] one? a positive error bound 6. 

Output: A cubic B-spline p(s) such that the approximate error between p(s) and r(t) is 

less than 6 and the approximate implicit spline for r{t). 

1. Divide the curve r(t) into quasi-cubic segments by Algorithm 3.15. 

2. Check the topology conditions. 

(a) Check the intersection of any pair of cubic Bezier curves which have the same 
Frenet frame at the endpoint, divide them to two parts on their shoulder points 
respectively, if they have common points more the endpoints. 

(b) Check the collision of any pair of tetrahedrons, divide them to two parts on their 
shoulder points respectively, if they have inner points. 

3. For each segment, find the cubic Bezier curves which approximate the given curve 
segment with precision 6 by Algorithm 4.3. 
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4. Find the implicit form for the cubic Bezier curves with the /i-basis method [31]. 

5. Convert the resulting rational cubic Bezier curves to a rational B-spline with a proper 
knot selection as the method presented in [2]. 

Remark 5.2. In the process of topology conditions checking, we only need to check the colli- 
sion of the sub-tetrahedrons subdivided from which are the intersected before the subdivision, 
since the sub-tetrahedrons are included in its father tetrahedrons. It means that the less and 
less pairs of tetrahedrons need to be checked in the subdivision process. 

Theorem 5.3. From Algorithm 5.1, we obtain a piecewise C^ continuous approximate cubic 
B-spline curve which keeps the singular points, inflection points, and torsion vanishing points 
of the approximated parametric curve. At cusps, the approximate curve is C° continuous. 

Proof. Algorithm 5.1 gives the G^ cubic Bezier spline since it is constructed as the hermite 
interpolation of the original curve, if the character points are not cusps. Then C^ continuity 
can be ensured from the conversion from the Bezier spline with a proper knot selection [2]. 
The singular points of the curve are treated as segmenting points. Since at the segmenting 
points, the left and right Frenet frames are preserved, the origin curve and the approximate 
curve have the same singular points. Since the cubic spline introduces no more singular 
points, the algorithm keeps the singular points. At a cusp, its left (right) tangent and 
osculating plane are kept according to Algorithm 3.15, and the approximate curve is then 
only C° continuous. 

The character points include the vertices of the topology graph. The topology conditions 
ensure that the topology is persevered while the topolgy line edges are replaced by the cubic 
Bezier curve segments. According to Theorem 4.6, the approximate curve from Algorithm 5.1 
converges to the approximated curve and they have the same topology. 

The left and right Frenet frames of the approximate curves are the same as that of 
the approximated curve at the character points, which means that the principal normal 
vector and the osculating plane are both kept. Then the principal normal vector changes its 
direction at the inflection point. Similarly, the curve does not pass through the osculating 
plane at the torsion vanishing point. D 

Finally, we give several examples to illustrate the algorithm. 

Example 5.4. The space curve ri(t) from Example 6 in [27] has a singular point (0,0,0) 
at t = ±1, ±oo, where 

ri(t) = 

The curve segment ri{t),t G [—2,2] and its approximate spline curve p(s) are shown in 
Figure 6, they are shown in the same figure for comparison and the tetrahedron sequence is 
also given in Figure 1, the numerical error e{t) is shown in Figure 8. 

As we know, the point (0, 0, 0) is a characteristic point from the topology determining. It 
is preserved in p(s) and p(s) is C^ at this point. Each corresponding segment of p{s) and 
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t'il-t') 
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'it' + lf' 


it' + I)' 





Figure 6: ri(i) and p(.s) 





Figure 7: ri(i) v.s. p(s) and eontrol tetrahedron 




Figure 8: Numerical error for ri with m — 300 

ri(t) is interpolated with the Frenet Frames at the endpoints. One can find that ri(t),t G 
[— oo, +oo] is an asymmetric space trijolium curve. To approximate the other two parts of 
t G [— C)0, —2] and t G [2, +oo], we can transform t = ±oo to t = by a reparametrization 
as t' = 1/t. Then approximating ri(t'),t' G [—1/2,1/2] and combining the former spline 
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segment, we can get the approximation of the whole trifolium curve. 

Example 5.5. Two more space curves are given in this example. r2(t) has a complex sin- 
gular point and r3(t) is a random curve with degree nine. 



r2(t) 



t^t 



t{t 



t{t 



t^f ' 1 + t2 



l + t2 



,tG [-1/16,3/2] 



r3(t) 



t(ll81 t8_i878 1''-1236 1''+1960 4^+2058 f*-2688 f^+532 1^-9+72 1) 
-2+9 1-72 (2+308 4^-840 t'i+1218 i^ -952 1^ +588 1^ -408 18+149 1^ 
t(-1686t^+287t*+3252t6_2464t^+462t"'+168t=' -28*2+9) 

-2+9 t-72 i2+308 t^-SAO t*+l2lS t5-952 16+588 1' -mS t« + 149 1^ ' 

4*2 (^263 4^-924 46+1338 15-1190 1"'+861 13_483 42+154 i-l8) 
-2+9 1-72 12+308 13-840 14+1218 15-952 46+588 i7_408 18+149 1^ 



,tG[0,l] 



The approximated curves, approximate spline curves, and the numerical errors are shown 
in the following figures (Figures 9, 10, 11). In r2{t), (0,0,0) is a self-intersected point with 
t = 0, 1, it is also a cusp point att = l. This point is preserved in our approximate B-spline 
curve p(s). Furthermore, the limited tangent directions of the cusp are also preserved. p(s) 
is C^ or C^ at (0, 0, 0) when p(s) passes through (0, 0, 0) as a self-intersected or a cusp point 
respectively. The approximation information for curves ri,r2, and r^ is listed in Table 1. 





Figure 9: Y2{t) v.s. P2(s) and r^lt) v.s. P3(s) 



curve 


degree 


error 


segments 


interval 


ri 


8 


0.004157 


8 


[-2,2] 


r2 


5 


0.0001677 


4 


1 3 
16' 2 


r3 


9 


0.03298 


6 


0,1] 



Table 1: Numerical Approximation 
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Figure 10: Numerical error for r2 




Figure 11: Numerical error for r^ 

6. Conclusion and further work 

We present an algorithm to construct a rational cubic B-spline approximation for a space 
parametric curve. The main purpose of the work is to present an isotopic approximation 
method which preserves the geometric features of the original curve. The approximated 
curve is divided into quasi-cubic segments which have similar properties to those of a cubic 
Bezier curve. Sufficient conditions are proposed for a divided segment having the expected 
properties and then its approximate Bezier spline is naturally constructed. Based on these 
properties, the shoulder point approximate algorithm is presented and it is proved to be 
convergent. An approximate implicitization can be found by the yU-basis method. The 
method is applicable for any parametric space curve in theory, although the given conditions 
are more difficult to compute when the parametric expression is not in rational form. 

The intersection curve of a parametric surface and an implicit surface is another impor- 
tant type of space curves. The curve can be regarded as parametric form with two parameters 
and a constraint function for them. As a further work, we will study the approximation of 
this type of space curve. 
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